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ABSTRACT 


The  calculated  stresses  and  displacements  induced  in  anisotropic 
plates  by  short  duration  impact  forces  are  presented  in  this  report. 

The  theoretical  model  attempts  to  model  the  response  of  fiber  composite 
turbine  fan  blades  to  impact  by  foreign  objects  such  as  stones  and  hail¬ 
stones.  In  this  model  the  determination  of  the  impact  force  uses  the 
Hertz  impact  theory.  The  plate  response  treats  the  laminated  blade  as 
an  equivalent  anisotropic  material  using  a  form  of  Mindlin's  theory 
for  crystal  plates.  The  analysis  makes  use  of  a  computational  tool 
called  the  "fast  Fourier  transform".  Results  are  presented  in  the 
form  of  stress  contour  plots  in  the  plane  of  the  plate  for  various 
times  after  impact.  Examination  of  the  maximum  stresses  due  to  im¬ 
pact  versus  ply  layup  angle  reveals  that  the  ±15°  layup  angle  gives 
lower  flexural  stresses  than  0°,  ±30°  and  ±45°  cases,  for  55%  graphite 
fiber-epoxy  matrix  composite  plates. 
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SUMMARY 


This  report  summarizes  the  recent  progress  in  the  attempt  to  model 
foreign  object  impact  of  fiber  composite  fan  blades  by  small  objects 
such  as  stones  and  hardstones.  The  high  speed  impact  of  these  objects 
with  composite  materials  should  result  in  short  duration  impact  times 
(<  100  psec) .  In  the  present  model  the  composite  blade  is 
represented  by  an  anisotropic  plate  of  infinite  extent.  Thus  only  the 
mechanics  in  the  area  of  the  impact  point  are  studied.  The  effect  of 
edges  and  stress  wave  reflection  are  not  treated  in  this  report,  though 
work  is  proceeding  on  this  aspect  of  the  problem. 

The  combined  impact  contact  dynamics  and  plate  response  are  separated 
into  two  sub-problems.  The  impact  force  time  history  is  determined  from 
the  Hertz  contact  theory  for  a  half  space,  while  the  plate  response  is 
determined  from  an  assumed  impact  pressure  distribution  over  the  impact 
contact  area.  Linear  elastodynaraics  is  used,  neglecting  viscoelastic- 
plastic,  nonlinear  and  fracture  effects.  Thus  it  is  hoped  that  the  model 
will  predict,  prefracture  or  predamage  stresses. 

Using  an  approximate  theory  of  anisotropic  plates  due  to  Mindlin, 
five  waves  are  shown  to  make  up  the  main  part  of  the  motion.  Two  of  the 
waves  involve  inplane  displacements  while  the  other  three  waves  involve 
flexural  plate  motion.  The  analysis  makes  use  of  a  computational  tool 
called  the  fast  Fourier  transform.  Several  computer  programs  were 
developed  to  calculate  the  stress  levels  behind  the  wave  fronts  due  to 
a  specified  impact  force  distribution.  The  output  of  these  programs  are 


in  the  form  of  stress  contour  plots.  It  was  generally  found  that  the 
stress  levels  were  highest  in  directions  along  the  fibers.  Also  the 
maximum  stresses  appear  to  propagate  in  the  lowest  flexural  mode. 

These  waves  are  found  to  be  highly  dispersive  and  change  their  shape 
as  the  wave  propagates .  Examination  of  the  maximum  stresses  due  to 
impact  versus  ply  layup  angle,  reveals  that  the  ±15° layup  angle  gives 
lower  flexural  stresses  than  0°,  30°,  and  45°  cases.  The  flexural 

stresses  for  the  15°  case  are  34%  lower  than  those  for  the  45°  layup 
angle  case.  For  the  interlaminar  shear  stresses,  the  values  seem  to 
be  insensitive  to  layup  angle.  For  the  average  membrane  stresses  due 
to  inplane  motion,  the  values  immediately  after  impact  appear  to  be 
lower  for  the  lower  fiber  layup  angles.  The  flexural  stresses  are 
found  to  depend  on  the  ratio  of  impact  object  radius  to  plate  thickness. 

Continued  work  is  in  progress  on  the  edge  impact  problem  and  the 
effect  of  leading  edge  protection. 
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INTRODUCTION 


This  research  has  been  motivated  by  recent  attempts  to  study  the 
impact  resistance  of  fiber  composite  materials.  These  materials  are 
being  considered  for  application  to  jet  engine  fan  or  compressor 
blades  and  must  withstand  the  forces  of  impact  due  to  the  ingestion 
of  objects  such  as  birds,  stones,  or  hailstones  at  speeds  up  to 
500  meters  per  second.  Recognizing,  of  course,  that  inelastic 
deformation  will  occur  at  the  impact  point,  we  none  the  less  pursue 
an  elastic  analysis  as  a  prelude  to  the  more  difficult  task  of 
predicting  permanent  damage.  In  this  report  we  are  interested 
principally  in  how  the  energy  propagates  awav  from  the  impacted 
area.  It  has  been  observed  that  damage  in  these  fiber  composite  plates 
occurs  away  from  the  impact  area  near  edges  and  boundaries  as  well  as 
at  the  impact  point. 

In  a  previous  paper  (ref.  1) ,  the  author  presented  a  mathematical 
model  for  stress  wave  propagation  in  anisotropic  plates  based  on  the 
work  of  Mindlin  and  co-workers.  Five  partial  differential  equations 
of  motion  were  obtained  for  orthotropic  symmetry  in  which  the  inplane 
and  flexural  motion  were  described.  The  two-dimensional  velocity  and 
wave  surfaces  were  presented  and  the  principal  vibratory  direction 
of  particle  motion  for  each  wave  normal  was  presented. 

Section  II  will  present  an  analysis  of  the  two-dimensional 
impact  problem.  This  analysis  is  based  on  the  use  of  a  Laplace 
•transform  on  time  and  a  two-dimensional  Fourier  transform  on  the 
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space  variables.  The  solution  permits  the  analytical  inversion  of 
the  Laplace  transform  while  a  computational  tool  called  the  "Fast 
Fourier  Transform"  (ref.  2)  is  used  to  numerically  invert  the  Fourier 
transform  solution.  Estimates  of  the  impact  time  and  force  are  made 
using  the  Hertz  impact  theory. 

As  a  preliminary  to  the  two-dimensional  problem,  the  author  tried 
this  analytical-computational  technique  on  a  few  simpler,  one -dimensional 
wave  problems.  This  section  will  present  the  results  of  that  study. 
Displacement  and  stresses  are  calculated  for  a  short  duration  line 
force  normal  to  an  anisotropic  plate.  The  responses  for  various 
fiber  layup  angles  are  compared. 
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SECTION  I:  ONE -DIMENSIONAL  TRANSIENT  WAVES  IN  ANISOTROPIC  PLATES 


The  mathematical  model  used  in  this  paper  is  called  an  effective 
modulus  theory  of  composites.  The  model  is  valid  so  long  as  the  scale 
of  the  changes  in  stress  levels  (e.g.  wavelength)  is  much  larger  than 
the  sizes  of  the  composite  constituents  (e.g.  fiber  diameter,  ply 
spacing,  etc.).  This  assumption  has  been  used  by  many  researchers 
to  derive  equivalent  elastic  constants  from  long  wavelength  or  long 
pulse  length  wave  propagation  tests  (e.g.  Tauchert  and  Moon,  (ref.  3)). 
Equivalent  moduli  determined  in  this  way  have  checked  very  closely  with 
elastic  constants  obtained  from  static  and  vibratory  tests.  Tauchert 
and  Guzelsu  (ref.  4) ,  have  used  ultrasonic  tests  to  measure  dispersion 
in  composites  and  found  no  significant  departure  from  the  effective 
modulus  theory  for  longitudinal  waves  in  boron  fiber/epoxy  matrix 
composites  until  the  frequency  exceeded  5  megahertz.  Shear  waves, 
however,  began  to  show  dispersion  at  about  0.5  megahertz.  This 
corresponds  to  a  wavelength  to  fiber  diameter  ratio  of  about  40. 

Sun  et  al .  (ref.  5),  in  calculation  on  laminated  composites,  reports 
similar  results.  Thus  the  effective  modulus  theory  used  in  this  paper 
is  subject  to  the  restriction  that  the  impact  pulse  spectrum  have  its 
most  significant  wavelengths  greater  than  about  100  times  the  fiber 
diameter.  Further  discussion  of  dispersion  in  composites  may  be  found 
in  a  review  of  the  subject  by  the  author,  (ref.  6) . 

Equations  of  Motion 

In  the  approximate  theory  of  elastic  plates  in  this  study  (ref.  7) , 
the  displacement  is  expanded  in  the  thickness  variable  by  using  Legendre 
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polynomials 


P„(n), 


u  =  P  (n)u^  (x  ,x  ,t) 

A.  -n  —  -n  ^  r\.  '>  ^  ^  ■' 


n=o  n 


(i-i) 


where  n  =  x  /b,b  is  the  half  thickness,  (see  Fig.  I-l)  .  A  variational 
2 

method  is  used  to  obtain  the  equations  of  motion.  For  the  purpose  of 
this  work,  the  series  was  truncated  at  n  =  2,  with  only  seven  terms 

..0  ..0  ..0  ..1  ..1  ,.(2)  TTa  nC*  oil  +*1^  a  e?  "f"  ■V'O  T  T 


used;  namely,  ^  .  This  leaves  all  the  strains, 

1  2  3  1  2  3  2 

except  e  ,  e  ,  as  linear  functions  of  the  thickness  variable  n 
12  32 

To  further  simplify  the  equations,  the  inertias  and  high  frequency 

r2') 

terms  (derivatives)  of  u^,  u^  were  dropped,  resulting  in  explicit 

2  2 

expressions  for  these  terms; 


9u0  9u0 

1  o  9x  o  o  9x  ^ 
1  3 


(1-2) 


q  -  2(C  1—  +  C 


12  9^  ^  ^23  9^  ^ 

I  3 


where  (-q  )  is  the  transverse  loading  on  the  top  plate  surface.  Elimination 
2 

r  2 1 

of  the  terms  u^  and  u^  ^  from  the  remaining  five  differential  equations 
2  2 

results  in  the  equations 


9^  .  92u§ 

9x2  ^  ^^55  ^13^  9x  9x  2C  ^ 

0^  I  3  22  I 


(I-3a) 


6 


8t' 


C 

33 


9x2 

3 


92u0 

+  C  - ^  +  (C  +  C  ) 

55  93^2  55  13 

1 


!!!L. 

9x  9x  ^  2C  9x 
1  3  22  3 


(I-3b) 


92u0 
_ 2 

9t2 


=  C 


66 


82u0 


^2uO 

_ ^ 

9x2 

3 


1 

"  b 


9x 


1 

"  b 


^^3 

9x 


y  + 


2b 


92uJ 

9t2 


2,,1 


=  c 


11 


9'^u 


9x^ 


92u1 


+  C 


1 


55  3x2 

3 


9"^3 

*■^55  ^  ^13^  9x3^, 


9t2 


-  ^  C 


^  92ui 

C  -  + 

3  3  9x2 

3 


66 


9u^ 


9x 


92ui 


55  3x2 
1 


3 

B 


-  ^  C 


f  3uO 


66 


9x 


u 


n 


■-12 


2C  9x 
22  1 

^  92uJ 

+  (C  +  C  )  .  .  ■ 

55  13 
55  15 

C32  9q2 

2C  W 
22  3 


u 


1  1 


where 


c 

=  C 

-  C2  /C 

11 

11 

12  22 

c 

=  c 

-  C2  /C 

33 

33 

32  22 

C 

=  C 

-  C  C  /C 

13 

13 

12  32 

C  =  kC 
66  66 


(I-4a) 


Cl-4b) 


(I-4c) 
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In  his  approximate  theory  of  crystal  plates  (ref,  7) ,  Mindlin  uses 
correction  constants  to  adjust  the  thickness  vibration  mode  to  match  the 

exact  theory.  Thus  he  replaces  C  and  C  by '  k  C  and  k  C  respectively. 

44  66  1  44  3  66 

In  our  case  C  =  C  and  k  =  k  =  Tr^/12  . 

44  66  13 

The  stresses  in  the  plate  for  orthotropic  symmetry  take  the  following 

form : 


11 


22 


33 


13 


8u 


3u° 

r  _  +  r  _ 

11  13  3x 

^  ^  b-^ 


3u°  .  3u^ 

■(  C  ^  +  C  3— 

33  13 


0 


3u^ 


3u 


1 


11  13  3x^ 


^  + 


X2 


9ul 

C  ^  +  C  ^ 
33  9x  13  3x 


9uj 


1 


=  C 


55 


3u} 

3u0 

^2 

'  "  - 

c 

f 

i  -  + 

\ 

9ui 

Jx~  * 

1  ' 

3x 

1 

55 

t  9x 

1  ' 

9x 

1 

^12  ""2 
22 


(1-5) 


C32  X2 

22 


The  interlaminar  shear  stresses  t  and  t  are  quadratic  functions 

12  32 

of  X  and  cannot  be  accurately  found  from  u  .  However  in  analogy 
to  classical  plate  theories,  we  can  integrate  the  original  stress 
equations  of  motion  to  obtain 


(1-6) 
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The  in-plane  stresses  t  ,  t  ,  t  are  comprised  of  two  parts. 

11  33  13 

One  part  is  uniform  through  the  thickness  and  is  governed  by  dispersion 
free  equations  (I-3a]  and  (I-3b).  The  flexural  part  of  the  stress  is 
linear  in  the  thickness  variable  and  is  governed  by  equations  (I-4a) 
through  (l-4c)  which  exhibit  strong  dispersion  of  wave  phenomena. 
One-Dimensional  Wave  Solutions 

We  imagine  an  infinite  plate  on  which  a  distributed  load  is  placed 
along  the  line  ji  •  ;^  =  constant  (see  Fig.  1).  The  vector  ]|i  is 
normal  to  the  line,  and  we  assumed  the  load  on  the  plate  is  independent 
of  the  position  along  the  line.  Thus  the  surface  load  has  the  form 

q(x  ,x  ,t)  =  q(]|i.:^,t) 

1  3  ^ 

Because  of  this  relation,  the  motion  is  also  assumed  to  have  this  form, 

i .  e . 


u°(x  ,x  ,t)  =  u°(n.r,t),  etc. 
113  1 


For  derivatives  of  functions  of  this  form,  the  following  expressions 
are  useful; 

— r -  =  f'cos  a,  — r -  =  f'sin  a 

9x  9x 

1  3 


9 


where  primes  indicate  differentiation  with  respect  to  the  variable 
C  =  n  •  r  ,  and 

ji-;^  =  X  cos  a  +  X  sin  a 

The  angle  a  ,  denotes  the  departure  of  from  the  symmetry 

axis  X  =  0  .  To  find  solutions  to  the  equations  (1-3)  and  (1-4), 

3 

one  writes  the  unknowns  in  the  forms 


uj  =  U^(ji.;r,t),  uO  =  U^(^:^,t) 


(1-7) 


uO  =  W(ji.^,t),  uj  =  -  b  'l'^(ji-^,t),  =  b  'l'^(;i.^,t) 


Note  that  q  is  a  prescribed  function  in  time,  along  a  line  parallel 
to  the  normal  n  .  This  reduces  the  differential  equations  of  motion 
to  equations  in  two  variables  ^  ,  t  where  ^  =  n  •  r  .  The  formal 
solution  to  this  initial  value  problem  is  found  by  taking  a  Laplace 
transform  on  the  time  variable  and  a  Fourier  transform  on  the  space 
variable.  These  operations  are  defined  as  follows: 


L[f]  =  f(s)  =  /  e  ^'^f(t)dt 
o 

CO 

T[f]  =  f(k)  =  /  e^^^f(Ocic 


(1-8) 


10 


Thus  when  the  displacements  have  the  forms  given  by  (1-7)  and 
the  operations  (1-8)  are  performed  on  the  equations  of  motion,  the 
following  algebraic  set  of  equations  result,  where  initial  values 
of  }i  are  assumed  to  be  zero  and  where  q  =  0  for  t  <  0  . 


— 

- 

— 

(A  +  ps^/k^) 

1 1 

A  cosasina 
12 

U 

1 

1 

-q  ika  cosa 

2  1 

k2 

A  coscxsina 

12 

(A  +ps2/k2) 
22 

0 

3 

-q  ika  sina 

2  2 

— 

— 

- 

- 

“  — 

(C  k2+ps2)  -ikC  cosa 

ikC  sina 

W 

a  /2b 

66  66 

66 

‘2 

pfl^ik  cosa  (A  k2+pQ2+pg2^ 

-k^A  cosasina 

T 

q  ka  cosa 

0  1 1  0 

12 

3 

2  3 

sina  -k^A  cosasina 

(A  k2+pSl2+pg2^ 

T 

-q  ika  sina 

0  12 

22  0 

1 

2  4 

, 

L  J 

-  — 

(1-9) 


(I-IO) 
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where 


A  =  C  cos^a  +  C  sin^a 
11  11  55 


A  =  C  sin^a  +  C  cos^a 
22  33  55 


A  =  C  +  C 
12  55  13 


=  3  C  /b2p 
0  66 


0  is  the  frequency  of  the  first  thickness  shear  mode  from  the  exact 


theory.  Also, 


a  =  C  /2C  ,  a  =  C  /2C  ,  a  =  C  /2bC  ,  a  =  C  /2bC 

1  12  22  2  23  22  3  12  22  4  23  22 


Inversion  of  Laplace  Transform 

When  qC?,t)  is  given,  these  equations  may  be  solved  in  the 
transformed  variables.  For  example  consider  the  following  case. 

CASE  I.)  Midplane  Motion;  Normal  Load  (q^  =  q^  =  0) 

0  =  iq  { (A  +s2/k2)a  -  A  a  sin^aicosa/kA 

1  2  22  1  12  2 

(I-ll) 

0  =  iq  { (A  +s2/k2)a  -  A  a  cos^alsina/kA 

3  2  11  2  12  1 


where 


A  =  (A  +  s2/k2) (A  +  s2/k2)  -  a2  cos^asin^a 
11  22  12 


These  solutions  in  the  transformed  plane  are  thus  of  the  form 


0  =  iq  P(s2]/kA(s^) 


12 


where  P(s2),  A(s2)  are  polynomials  in  .  Afs^)  has  four  zeros  in  the 

complex  s  plane  for  each  k.  The  two  roots  Is  /k|  and  |s  /k| 

1  2 

correspond  to  the  two  real  wave  speeds  associated  with  the  propagation 
of  plane  waves  in  the  plate  (ref.  1).  In  fact,  ,  is  simply 

the  two-dimensional  acoustic  tensor  for  the  plate  as  discussed  in  a 
previous  paper.  We  are  thus  guaranteed  four  pure  imaginary  roots 
for  s 

s  =  ±iv  k  ,  ±iv  k 

i  2 


The  inversion  can  then  be  done  by  use  of  the  convolution  theorem; 
i.e. ,  if 


G(t)  =  L 


P(s^)  j 
A(s2) 


then 


U(k,t3  =  L'^[6]  =  i  q  (k,t-T)G(T)dT 


(1-12) 


For  our  case  it  can  easily  be  shown  that 


G(t) 


kP(-k2vp  sinkv^t  kP(-k2v2)  sinkv^t 


V 


(v2-v2) 
1  2 


(v2-v2) 
1  2 


CASE  II.)  Flexural  Motion;  Normal  Load. 


(1-13) 


In  a  similar  fashion  when  q  is  given,  W,  ?  or  T  can  be 

1  3 

solved  and  have  the  form 
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(1-14) 


R(s^) 

A  (s^jk^) 

1 


However,  in  this  case  the  roots  of  A  =  0  are  not  proportional  to  k  . 
This  means  that  the  velocity  of  wave  propagation  depends  on  the  wavelength. 
As  k  it  is  known  that  the  wave  speeds  become  independent  of  k  and, 

in  fact,  two  of  them  equal  the  non-dispersive  wave  speeds  v^ ,  v^  found 
in  the  previous  example  (ref.  1) .  It  is  also  known  from  the  dispersion 
relations  for  plates  (ref.  7),  that  for  k  real,  there  will  be  six 
pure  imaginary  roots  of  A(s2)=0at 


s  =  ±ia)  (k)  ,  iio)  (k) ,  iiw  (k) 
1  2  3 


The  relations  oj^(k)  are  known  as  the  branches  of  the  dispersion  relation 
for  these  plates.  One  branch  goes  through  the  origin,  i.e. 

m  (k)  0  as  k  ^  0 

1 

The  other  two  branches  are  sometimes  called  optical  since  as  the 

wavelength  becomes  infinite,  i.e.  k  ^  0,  w  and  w  denote  the  vibra- 

2  3 

tory  frequency  of  the  thickness  shear  mode  for  the  plate 

oj  (0)  =  0)  (0)  =  Q 
2  3  0 

A  typical  dispersion  curve  is  shown  in  Figure  2. 

The  inversion  again  makes  use  of  the  convolution  theorem  and  the 

residue  calculus  to  invert  1/A  (s^)  .  Thus  we  obtain 

1 
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(1-15) 


W(k,t)  =  q(k ,t-T)H(T) dx 

0 

RCto^)  sinw^t 

12  3  2 

P,(to  )  sina),t 

_ 3 _  _ f_ 

13  2  3 

Numerical  Inversion  of  the  Fourier  Transform 

It  is  at  this  juncture  that  the  solutions  to  most  problems  in 
transient  wave  propagation  are  limited  by  analytic  skill  in  extracting 
information  from  the  inversion.  In  a  few  cases  the  complete  solution 
can  be  given,  while  in  most,  far  field,  near  field,  short  and  long 
time  approximations  must  be  invoked.  The  numerical  inversion  of  the 
transform  has  until  recent  years  been  limited  by  the  calculation  of  an 
integral  for  each  point  in  the  space  ^  .  However,  recently  a  technique 
has  been  developed  which  obviates  the  need  for  a  separate  inversion 
for  each  point  in  the  real  space  t,  .  Known  as  the  "Fast  Fourier 
Transform"  (ref.  8) ,  it  takes  a  sampling  matrix  of  the  transform  of  U, 
say  (U(k  ),  U(k  ),  ...  U(kj^j))  and  returns  a  sampling  matrix  of  the 
original  function 

(U(?  ),  U(?  ),  U(c  ),  ...  U(c  )) 

12  3  ^ 

For  a  large  enough  N,  one  will  obtain  a  picture  of  the  original  function 
U(c).  Application  of  this  algorithm  to  our  problem  is  described  below. 


H(t)  = 


R(a)j) 


((Jj2-(jj^)  ((JJ^-OJ^) 

2  1  3  1 


sincjjt 

0) 
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The  specific  operation  that  these  computational  algorithms  perform 
is  called  a  finite,  discrete  Fourier  transform.  Thus  if  Dfl)  is  a 
one  X  N  matrix  of  data,  the  output  of  the  algorithm  T(J)  is  given  by 


T(J) 


(1-16) 


The  properties  of  this  operator  are  similar  to  the  conventional  Fourier 
operators  in  the  space  of  continuous  variables  (ref.  8) .  The  sum  is 
not  performed  in  a  direct  manner  on  the  computer  but  in  a  manner  which 
reduces  the  number  of  arithmetic  steps  and  makes  the  operation  attractive 
as  a  computational  tool. 

In  the  solution  of  partial  differential  equations  by  Fourier 
transforms  we  are  led  to  evaluate  integrals  of  the  type 

U(?)  =  ^  r  U(k)e""’"^dk 

—  CO 

If  significant  changes  in  U(c)  occur  over  distances  greater  than 
A  ,  then  the  largest  wave  number  of  interest  will  be 

K  =  2-n/X 

Thus  we  may  be  satisfied  with  an  approximation  to  U(C)  of  the  form 

U(?)  =  ^  r  U(k)e~^^^dk 

-K 
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or  shifting  the  coordinates 


itc^ 

e _ 2 

2:7 


UCk-K)e'^’^^dk 


This  latter  integral  may  further  be  approximated  by  the  limit  of  a  sum 


U 


K  iKt 
ttN  ® 


E 


U(kj-K:)e 


(1- 


where 


or 


^  = 


2k 


Ci-1) 


4tt 

NX 


(I-l) 


Aj  l<- 

U  =  — 

ttN 


N 

E 

1=1 


U(kj-K) 


-1217(1-1)  2? 


N 


X 


So  far  ll(c)  has  remained  a  continuous  function  of  Z  .  However  at 
the  points 


=  |(J-1)  ,  J  =  1,  2,  ..  N 


the  summation  becomes  identical  to  the  finite  descrete  Fourier  transform 


T(J) 


(1-18) 


iJ(Sj)  =  U(J) 


iTr(J-l)  (N-1)  K 
®  N  ttN 


where  T(J)  is  defined  in  (1-16)  and 


D(l)  =  U(kj-K) 


Note  that  this  scheme  gives  details  of  U(c)  of  distance  no  smaller  than 
X/2  .  This  was  implicit  in  choosing  maximum  wave  number.  Further  when 
N  is  even,  the  output  matrix  U(J)  describes  the  function  only  up  to 
_  1)  .  When  ij(kj-K)  is  real,  and  N  even,  we  have  the 

identities 


?)  (|  +  1)  =  0  ,  lmag(U(J))  =  0 

N  N 

U(y  +  1  +  k)  =  -  U(2  +  1  -  k) 


Thus  the  one-dimensional  output  matrix  il(j)  which  approximates  our 
original  function  is  antisymmetric  about  J  =  j  +  1  and  symmetric 
about  J  =  1  .  Since  U(J)  =  0  for  J  =  1  +  N/2  we  must  choose 
N  large  enough,  or  the  time  small  enough  to  ensure  that  our  wave  has 
not  reached  <;  =  A(N  +  2)/4  .  In  this  approximation,  we  have 
effectively  replaced  a  single  impact  source  by 
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an  infinite  periodic  array  of  sources  and  negative  sources.  Our 
approximation  will  be  valid  as  long  as  the  waves  from  each  of  these 
sources  do  not  interact  with  each  other. 

Notes  on  the  Numerical  Fourier  Inversion 

There  were  several  checks  made  of  the  fast  Fourier  computer 
routine  (ref.  8)  used  in  this  paper.  First,  known  functions  were 
transformed  analytically  and  inverted  numerically.  The  results  of 
these  tests  revealed  that  the  inversion  program  works  best  on  continuous 
functions  if  one  is  to  avoid  spurious  oscillations  near  points  of 
discontinuity . 

Second,  the  initial  value  problem  for  a  nondispersive  medium 
such  as  a  string  was  programmed.  The  results  of  this  test  are  shown 
in  Figure  3a.  The  shape  of  the  string  at  the  initial  time  and  subsequent 
times  as  predicted  by  the  computer  preserve  the  original  triangular 
shape  and  exhibit  the  wave  shift  at  the  proper  speed. 

Lastly  the  impact  of  a  string  was  programmed  using  the  same  force 
applied  to  the  plate  problem.  The  theoretical  result  predicts  a  constant 
displacement  behind  the  wave  front  as  shown  in  Figure  3b. 

For  the  line  impact  of  an  anisotropic  Mindlin  plate,  a  specific 
load  distribution  was  chosen  for  ease  of  analytical  calculation  of  the 
Fourier  transform  and  Laplace  convolution  integral; 


for 


for 


X,  <  a  and  t  <  t 


n  1 

%  '  ft'  ^ 

o 


o  2 

^  I  >  a  or  t  >  T 


(1-19) 


0 
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By  using  the  fast  Fourier  transform  (with  a  proper  sign  change) 
one  could  find  the  transform  for  an  arbitrary  load  distribution. 

However,  this  has  not  been  done.  Nor  has  the  author  used  the  load 
distribution  for  a  Hertz  contact  problem,  as  might  be  supposed  in  an 
impact  problem  (reasons  for  this  will  be  given  below) .  However  the 
chosen  distribution,  it  is  hoped,  will  exhibit  most  of  the  salient 
effects  of  the  impact  of  a  plate. 

While  ad  hoc,  the  particular  choice  of  the  load  distribution  is 
not  entirely  arbitrary.  Continuity  of  the  first  and  second  derivatives 
(as  the  function  (1-19)  exhibits)  is  dictated  by  a  desire  to  have  all 
stresses  continuous  (i.e.  to  avoid  shocks)  and  thus  snurious  oscillations 
in  the  numerical  inversion.  This  stems  from  the  fact  that  in  the 
Mindlin  theory  the  midplane  stresses,  having  as  wave  sources  the 
term  Vq2  ,  have  the  following  form  for  their  transforms 

^  rt 

t  ^  J  kq  (k,T)  sinkv(t-T)dT  ,  (a,B  =  1,3) 

“3  o  2 

where  v  denotes  either  the  first  or  second  wave  speed.  Thus  when 
has  the  form 

q  =  Q(k)6(T) 

2 


the  stresses  are  proportional  to 


0,  k  Q(k) 


sinkvt  V  T 


[F(vt-0] 


wherOj 


F(vt-?)  =  1,  <  vt 

=0  1^1  <  vt 
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Thus  for  continuous  stresses  at  early  times  after  the  wave  arrival,  the 
spacial  part  of  must  have  continuous  second  derivatives,  which  led 

to  the  choice  of  (1-19)  .  This  conclusion  was  found  also  in  the  numerical 
results  when  non-smooth  load  distributions  were  used. 

Results  for  the  Line  Impact  Problem 

Using  the  transient  load  distribution  described  above,  calculations 
were  made,  on  an  IBM  360-91  computer,  of  the  stresses  and  dis¬ 
placements  in  both  a  classical  and  Mindlin  plate  due  to  impact  type 
loading.  The  results  are  shown  in  Figures  4-8. 

One  of  the  important  parameters  in  this  problem  is  the  ratio  of 
load  extent  to  plate  half  thickness,  a/b.  When  a/b  is  of  order 
unity  or  less,  one  would  expect  that  shear  effects  would  become  important 
and  the  Mindlin  and  classical  plate  solutions  would  differ.  This  is 
so  ,as  shown  in  Figure  4  for  the  case  of  a/b  =  1  .  However  in  Figure  5 
when  a/b  =  10  the  displacements  calculated  by  both  Mindlin  and  classical 
theories  do  not  differ  by  very  much. 

The  remaining  figures  are  for  the  Mindlin  plate  and  were  calculated 
for  the  composite,  55%  graphite  fiber-epoxy  matrix,  using  equivalent 
elastic  constants  for  various  layup  angles  of  the  fibers.  The  constants 
were  taken  from  Chamis  (ref.  9) . 

In  Figure  6  we  have  plotted  the  center  plate  displacement  versus 
time  for  both  the  classical  and  Mindlin  theories.  For  the  classical 
plate  this  function  can  be  found  explicitly.  The  displacement  increases 
as  the  square  root  of  time  when  the  impact  force  is  a  delta  function 
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in  time  and  space.  This  result  is  confirmed  by  the  numerical  results 
in  Figure  6  and  is  also  the  case  for  the  Mindlin  plate  for  large  time. 
(Note  that  for  a  string,  the  displacement  at  the  center  is  time 
independent  after  impact.) 

The  effect  of  layup  angle  on  the  plate  response  is  shown  in 

Figure  7,  for  various  times  after  impact.  As  can  be  seen,  the 

displacement  is  somewhat  insensitive  to  layup  angles  of  up  to 

about  ±15°  for  either  line  loads  along  the  x  axis  or  along  the 

3 

X  axis . 

1 

Finally  in  Figure  8  the  induced  membrane  or  average  in-plane 

stress  1/2 (t  +  t  )  is  shown  for  various  layup  angles  at  a 
11  33 

time  immediately  after  the  impact  time.  One  can  see  that  the 
initial  compression  pulses  are  proceeded  by  wave  fronts  which 
vary  in  magnitude  and  velocity  with  layup  angle.  Also  for  a 
wave  propagating  along  the  x  axis  (load  on  the  x  axis) , 
the  ±45°  layup  results  in  higher  stresses  than  for  other  layup 
angles.  While  the  flexural  stresses  are  higher  than  the  membrane 
stresses,  the  membrane  stresses  will  propagate  ahead  of  the  bending 
waves  and  will  have  a  tension  pulse  in  the  signal  which  might 
cause  splitting  through  the  plate  or  plys  as  have  been  observed 
in  experiments . 

In  the  case  of  one -dimensional  wave  propagation  in  a  direction 
n  =  [cosa,  sina],,  the  direction  of  the  displacement  is  not  parallel 
to  but  is  known  for  the  particular  wave  in  question.  The  direction 
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of  displacement  for  the  first  two  in-plane  waves  was  given  in  Figure  6 
of  Ref.  1 


=  UCjT.‘;^-vt)  [cos  g,  sin  6] 


From  (1-5),  the  stress  may  be  determined.  In  particular,  the  average 

in-plane  stresses  t^  ,  t°  ,  t°  can  be  found  as  functions  of  the 

11  33  13 

mean  stress  a  =  1/2 (t°  +  t°  ) ; 

1 1  33 


t°  =  2a(C  cos  a  cos  g  +  C  sin  a  sin  g)/D 
11  11  13 


t*^  =  2a (C  sin  a  sin  g  +  C  cos  a  cos  g)/D 

33  33  13 


t*^  =  2  a  C  (cos  g  sin  a  +  cos  a  sin  g)/D 

1 3  55 


where 


D  =  (C  +C  )  cos  a  cos  g  +  (C  +  C  )  sin  a  sin  g 
11  13  33  13 
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SECTION  II:  STRESS  WAVES  DUE  TO  CENTRAL  IMPACT  OF  ANISOTROPIC  PLATES 


In  Section  I  of  this  report,  an  analytical-computational 
method  was  described  for  the  calculation  of  the  impact  induced  stresses 
behind  the  wave  front.  The  method  was  checked  for  known,  one-dimensional 
problems  such  as  the  impact  of  a  string,  and  the  line  impact  of  a  classical 

anistotropic  plate.  The  induced  stresses  were  given  in  terms 

of  the  impact  pressure.  However,  no  prediction  was  made  of  the  maximum 

impact  pressure  in  terms  of  the  velocity  and  mass  of  the  impacting  body. 

In  this  part  of  the  report,  the  problem  of  normal  or  central  impact 
is  discussed.  Two-dimensional  stress  patterns  are  presented  in  terms  of 
the  impact  pressure  and  an  estimate  is  made  of  its  magnitude  for  a 
spherical  impacting  body  of  known  velocity  and  mass. 

Description  of  the  Mathematical  Model 

Before  discussing  the  results,  a  discussion  of  the  assumptions  made 
in  the  model  will  be  given.  In  practice,  fiber  composite  plates  are  made 
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up  of  a  number  of  unidirectional  plys  oriented  at  various  angles  to 
obtain  certain  desired  properties.  When  the  properties  and  lay  up 
angles  of  each  ply  are  known,  one  can  obtain  the  plate  constants  by 
averaging  across  the  thickness.  These  average  elastic  constants 
represent  an  equivalent  anisotropic-homogeneous  material.  In  the  present 
work  the  laminated  plate  is  replaced  by  an  equivalent  anisotropic  plate. 

The  equivalent  elastic  constants  were  obtained  from  a  computer  code 
developed  by  Chamis  (ref.  9) .  The  following  additional  assumptions  are 
implicit  in  this  model. 

Linear  Elastic  Properties 

This  model  is  based  on  linear  elasticity.  Thus,  plastic,  viscoelastic, 

and  fracture  effects  are  not  taken  into  account.  The  results  of  these 

effects  is  to  decrease  the  amount  of  wave  energy  that  can  propagate  into 

the  plate.  Thus,  the  elastic  analysis  represents  an  upper  bound 

on  the  actual  stresses  during  impact.  The  effect  of  initial  stress  in 
the  plate  has  also  not  been  considered  in  this  report. 

Boundaries  and  Tapering 

Another  limitation  of  this  work  is  the  neglect  of  the  effect  of 
boundaries  or  edges  of  the  composite  fan  blade.  The  reflections  from 
boundaries  can  of  course  be  handled  with  the  linear  elastic  theory. 

However,  if  one  is  interested  in  the  maximum  stresses,  these  should  occur 
near  or  at  the  impact  point  except  for  the  case  of  edge  or  near  edge 
impact.  Also,  for  small  impact  times,  the  waves  behave  as  if  the  finite 
plate  were  infinite.  The  case  of  edge  impact  is  under  study  at  this 
writing . 

In  this  report  the  plate  thickness  is  assumed  to  be  constant,  while 


25 


in  a  fan  blade  taper  and  twist  will  be  present.  The  neglect  of 
these  effects  seem  to  be  of  lesser  importance  than  those  due  to 
inelasticity  and  boundaries. 

Thickness  Reflections  and  the  Mindlln  Theory 

The  use  of  an  equivalent  anisotropic  plate  neglects  the  reflections 
of  waves  at  the  ply  boundaries.  However,  this  approximation  will  be  valid 
for  wavelengths*  greater  than  the  ply  thickness,  and  valid  also  for 
wavelengths  greater  than  the  plate  thickness.  There  are  analyses  which 
consider  wave-ply  Interactions,  (Ref.  10),  but  in  general  the  change 
in  stress  must  occur  over  distances  comparable  to  ply  or  fiber  thickness 
for  these  effects  to  become  important. 

In  the  Mindlin  theory  of  plates,  the  wave  reflections  are  averaged 
through  the  thickness,  (ref.  7).  In  dynamic  loading  of  a  plate,  the  stresses 
will  propagate  through  the  plate  thickness  as  well  as  away  from  the 
impact  point.  These  waves  will  suffer  many  reflections  as  they  propa¬ 
gate  back  and  forth  between  the  plate  surfaces.  The  calculations  of  these 
many  reflections  for  a  three-dimensional  plate  thus  become  impos¬ 
sible  for  anything  but  very  early  times  after  impact.  The  Mindlin  plate 
theory  thus  restricts  the  mathematical  equations  to  a  description  of 
the  average  plate  midplane  motion  and  rotations. 

In  the  Mindlin  theory  used  in  this  report  the  plate  displacements 
u  =  (uj^,  u^,)  (See  Fig.  9)  have  the  form  (See  Ref.  1) 


*An  effective  wavelength  can  be  defined  as  the  contact  time  of  impact 
multiplied  by  the  smallest  wave  velocity  in  the  material. 
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Cii-1) 


The  equations  governing  the  motion  are  listed  below.  The  first  two 
equations  are  for  the  inplane  displacements  u°,u°  and  the  next  three 
govern  the  transverse  displacement  u°  and  flexural  rotations  uj^/b  , 
u^/b  .  The  transverse  impact  loading  is  represented  by  -q„(x, ,x„,t) 
which  we  will  prescribe  in  the  following  section. 


9^uJ  9^uJ  9^u®  8^u^ 

P - C  - +  c  — -  +  (C  +  C  )  — ^ 

^  ,-2  11  :i„2  55  .,,2  55  13  „ 


"l2  '^2 


3xi  8x3  2C22  8X3 


(11-2) 


82u0  82^0  82^0 

p  -  =  C33  - ^  +  C  - 

8  t^  8x2  8x2 

3  1 


82u0  82u0  C  Sq. 

— Z  +  (C„  +  C.3)  - ^ - +^-i- 

8x2  8X3  8x2  ^^22  ^^3 


82u'^  82u°  .  8uj  1  ^^3  1 

c. .  — ^  +  C,  ,  - ^  +  C, ,  ^  ^  +  C .  .  ^  qo 

66  3^2  44  3^2  66  b  44  b  3^  2b  2 

1  3  1  3 


8x2 


82u^  82u3 

-7  ^  <=55  — 


8x  8x 
1  3 


ic 

^  ^  ^S2 


(11-3) 


+  C„ - 1  +  (C33  +  C33) 


33  ,9  ^55  ,  0 


8x  8x 
1  ^ 


o  8u®  u^  C„„  8q 

_  3  „  ,  2  ,  3  >,  +  23  2 

b  ^44''  -  3-  —  )  + 

8x  b  2C„„  8xp 

3  22  3 
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where 


Si  =  Si 


C  =  C 
33  33 


12 


22 


32 


22 


C  C 

r  -  r  -  C  =  kC 

13  13  n  ’ 

^22 


The  factor  k  is  a  correction  factor  introduced  by  Mindlin  and 
equal  to  k  =  Tr^/12 

If  the  plate  is  struck  at  a  point,  the  energy  will  propagate  into 
the  composite  in  the  form  of  elastic  waves  (see  Ref.  1)  .  Two  of  these 
waves  will  have  anisotropic  wave  fronts  and  will  depend  on  the  fiber 
layup  angle  of  the  composite  as  shown  in  Figure  10.  The  average  in-plane 
stresses  across  the  thickness  will  propagate  in  this  manner.  The 
flexural  energy  will  nropagate  behind  a  slower  isotropic  wave  front. 

The  stresses  in  the  plate  are  given  in  terms  of  the  displacement 
(Section  I  of  this  report) .  The  mathematical  problem  to  be  solved 
consists  of  the  following:  find  a  solution  to  the  coupled  partial 
differential  equations  (II-2),  (II-3)  when  the  loading  function  is  a 

prescribed  function  of  the  plate  coordinates  and  time,  and  the 

plate  is  initially  at  rest. 

The  solution  of  this  problem  was  accomplished  using  a  combination 
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analytical  and  computational  technique  involving  Laplace  and  Fourier 
transforms.  Details  of  finding  the  solution  are  discussed  in  Appendix  A 
of  this  report.  Results  of  these  calculations  are  discussed  below. 

Before  considering  the  results,  we  must  discuss  the  choice  of  the 
impact  loading  function  q2  • 

To  determine  the  maximum  contact  pressure  one  must  know  the 
impact  pressure  distribution  on  the  plate  surface.  In  the  static  Hertz 
theory,  this  pressure  is  given  by  (see  Ref.  11  and  12) 

q  =  -P  (1  -  (-)2)1/2  (II-4) 

^2  0  a 

This  distribution  is  not  suitable  for  our  theory  since  the  infinite 
slope  of  P(r),  at  r  =  a  ,  would  appear  in  our  approximate  equations 
(II-2)  and  (II-3)  as  a  source  term.  In  fact,  it  is  shown  in  Section  I  of  this 
report  that  the  second  derivative  of  P(r)  must  be  finite  at  r  =  a 
for  the  stresses  to  be  continuous.  To  resolve  this  problem,  an  ad  hoc 
pressure  distribution  is  assumed  which  produces  finite  stresses; 


q  =  -P  (1  -  2(^)2  +  r  <  a 

2  0  a  a  — 

q  =  0  ,  r  >  a 
2 


(II-5) 


f(t)  =  sin  Trt/T  ,  t  1,  T  ,  f(t)  =0  t  >  i 
It  should  be  noted  that  in  an  actual  impact,  the  impact  area 
changes  with  time.  This  effect  is  difficult  to  model  in  the  analytical 
part  of  the  problem.  Instead  we  have  chosen  to  prescribe  a  time  varying 
pressure  over  the  maximum  contact  area  . 
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The  total  force  produced  by  such  a  pressure  distribution  is 


r  3- 

F  =  2Tr\  q(r)  r  dr 

°  'o 

ra^P 

TT  -  _ Q- 


Uslng  (  9  )  we  may  calculate  the  total  impulse  to  be 

-  T 

I  =  i  F  sin—  dt  =  4  a2p  t 
0  T  3  0 

'  0 


CII-6) 


(11-7) 


If  the  rebound  velocity  is  V  ,  and  the  inital  velocity  V  , 

1  0 

the  impulse  given  by 

I  =  M(V  +v  )  (11-8) 

1  0 

Herztian  Impact  Stress 

In  order  to  determine  the  impact  stress  equations  (II-2,  3),  the 
impact  pressure  distribution  must  be  found  in  terms  of  the  object  mass, 
velocity,  geometry  and  elastic  properties.  No  exact  solutions  to  this 
dynamic  problem  are  known  in  the  theory  of  elasticity.  However  the 
static  theory  of  contact  of  elastic  isotropic  bodies  is  known  and  was 
developed  by  Hertz  (ref.  11).  It  is  known,  however,  that  even  under 
small  contact  forces  the  contact  pressure  exceeds  the  elastic  limit 
of  conventional  materials  (see  e.g.  Goldsmith,  (ref.  12)).  One  would 
expect  then  that  for  high  speed  impact  (>100  m/sec.)  the  local  contact 
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stresses  would  exceed  the  elastic  limit  of  most  materials.  In  suite  of 


this  violation  of  the  elastic  assumption  of  the  Hertz  theory  of 
impact,  experiments  have  shown  that  the  impact  times  predicted 
by  the  elastic  theory  agree  reasonably  well  with  data  from  impact 
tests.  The  stresses  predicted  by  this  theory  however  are  upper 
bounds  on  the  actual  stresses.  Also  if  the  surface  of  impact  belongs 
to  a  structure  which  can  move,  then  the  contact  times  predicted  by 
this  theory  are  lower  bounds  on  the  actual  time  of  contact  of  object 
and  structural  surface.  With  these  restrictions  in  mind  the  Hertz 
theory  of  impact  will  be  reviewed  as  it  can  be  adapted  to  composite 
materials . 

Consider  the  contact  of  a  sphere  of  radius  R  and  elastic  half 
space.  The  contact  force  between  the  two  bodies  F,  is  related  to 
the  relative  approach  of  the  bodies  a,  and  has  the  dependence 


F 


(11-9) 


where  k  is  a  constant  dependent  on  the  properties  of  the  bodies. 
2 

When  both  bodies  are  isotropic  this  constant  is  given  by 


k 

2 


4  .1/2 
3 


l-v2  l-v2 

^  "E 

1  2 


-1 


(II-IO) 


where  v.  Fare  Poisson's  ratio  and  Young's  modulus  respectively  for 


each  body. 

Composite  materials  however  are  anisotropic  in  general.  The 
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corresponding  problem  for  the  contact  of  anisotropic  bodies  has 

recently  been  given  by  Willis  (ref.  13) .  For  this  case  the  contact 

area  is  an  ellipse,  the  dimensions  of  which  must  be  determined  from 

algebraic  equations .  The  shape  of  this  ellipse  for  typical  composite 

anisotropy  has  not  been  determined,  though  from  experiments  of  the 

author  on  graphite/ epoxy  and  boron/aluminum  composites  seem  to  show 

that  the  contact  area  is  close  to  a  circle.  It  is  assumed  then  that 

for  the  impact  of  an  isotropic  sphere  with  a  composite  surface  the 

contact  area  is  a  circle.  Also  the  constant  k  is  replaced  by 

2 


k 

2 


l-vf 


(11-11) 


where  E  ,v  are  the  elastic  constants  of  the  sphere,  and  C  is  the 
11  22 

elastic  modulus  of  the  composite  plate.  This  assumption  of  course  is 
just  an  educated  guess. 

In  the  Hertz  impact  theory  the  force  (11-9),  which  was  determined 
from  a  static  solution,  is  used  in  Newtons  law  for  the  sphere 


3/2 

k  a 
2 


(11-12) 


where  M  is  the  mass  of  the  sphere  and  V  the  instantaneous  velocity  of 
the  sphere . 

This  assumption  is  only  valid  when  the  contact  time  is  much  greater 
than  the  time  for  elastic  waves  in  the  sphere  to  traverse  the  object. 
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Further  when  the  motion  of  the  plate  becomes  large  during  the  contact 
interval,  (11-12)  must  be  replaced  by 


=  k(U-W)^/^  (11-13) 

dt2 


where  U  is  the  displacement  of  the  sphere  and  W  is  the  displacement 
of  the  plate  at  the  point  of  contact.  This  problem  requires  simultaneous 
determination  of  the  motion  of  the  sphere  and  plate.  This  has  not  been 
done  in  this  report  but  work  on  this  problem  is  in  progress. 

IVhen  W  Ih  0  during  the  contact  period,  the  time  of  contact  can  be 
shown  to  be  (see  e.g..  Goldsmith  (ref.  12)) 


T 


da 


[V2  - 


M 


,,5/2]  1/2 


2.943 

_ 


(11-14) 


where  a  is  the  maximum  approach  of  the  impacting  object.  This  value 
is  given  by 


a 

1 


5 

4 


2/5 


MV2 

k 


2 


(11-15) 


For  elastic  impact  the  maximum  radius  of  contact  is  given  by 


a 


(11-16) 


The  impact  force  in  this  theory  as  determined  from  the  solution  of 
(11-12)  is  given  by 
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,  t  <  T 


F  =  F  sin  TTt/T 
0 


=  0 


,  t  >  T 


(11-17) 


where  F  =  3.36  —  .  The  pressure  distribution  in  this  theory  is  given 

0  T 

by  tlie  expression 

r  .  1/2 

P  =  P  (1  -  (-)2) 

0  ^ 

'7 

where  P  =  -pit  F  .  It  is  interesting  to  note  that  the  impact 

0  2  0 

pressure  P  ,  is  independent  of  the  radius  of  the  impacting  object. 

0  - - 

This  distribution,  however,  is  determined  from  a  static  elasticity 
solution  and  should  not  be  expected  to  hold  under  dynamic ■ impact .  Of 
course  the  contact  radius  varies  with  time  reaching  a  maximum  value 
(II-4)  at  half  the  period. 

Another  limitation  of  the  theory  is  the  fact  that  it  predicts 
a  rebound  velocity  equal  to  the  initial  velocity  .  In  an  actual 
impact,  momentum  would  be  transmitted  to  the  structure,  thus  changing 
rebound  velocity  to  something  less  than  V  .  Solution  of  the  coupled 
problem  (11-13)  should  enable  a  better  determinatinn  of  the  reboimd 
velocity. 

A  further  improvement  of  the  theory  might  be  achieved  if  a  more 
general  contact  law  is  used,  say  (called  the  Meyer  Law,  see  (ref.  12)) 

F  =  k  a 
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where  n,  k,  would  be  determined  from  static  experiments  on  comnosite 
materials.  This  approach  is  presently  under  investigation. 

Well  aware  of  the  long  list  of  limitations  on  the  above  model  for 
impact,  we  have  nonetheless  used  the  above  procedure  to  obtain  estimates 
of  the  impact  time,  contact  force,  and  maximum  contact  pressure  for 
different  impacting  bodies  and  speeds.  The  results  of  these  calculations 
are  shown  in  Figures  11  and  12  where  data  is  presented  for  spherical 
ice  particles  and  for  granite-like  stones. 

For  impact  speeds  from  100  to  500  meters/sec.  and  1/2  to  2  cm. 
diameter  granite  spheres,  the  contact  times  range  from  15  to  85  ysec. 

In  summary,  these  impact  formuli  reveal  the  following  dependence 
on  impact  velocity: 


T  -v  1/V 


1/5 


P  'v  V 


2/5 


These  results  should  only  be  used  as  guidelines,  since  many  assumptions 

are  made  which  break  down  at  high  velocities. 

Goldsmith  and  Lyman  (ref.  14)  have  shown  the  Hertz 

theory  to  be  remarkably  valid  insofar  as  contact  time  and  peak  pressure 

are  concerned  for  the  impact  of  hard  steel  spheres  (1/2”  diameter)  into 

a  hard  steel  surface  for  velocities  up  to  300  ft. /sec.  (91.5  meters/sec.). 

The  calculations  in  this  report  based  on  the  Hertz  law  of  impact  extend 

well  beyond  this  limit  (100-500  meters/sec.).  Thus  the  contact  times 

and  maximum  impact  pressures  presented  in  Figures  11>  12  for  graphite/epoxy 

can  only  be  used  as  rough  guides  for  which  the  values  for  contact  times 

are  lower  bounds  on  the  actual  times  and  the  values  for  pressure  are 

upper  bounds . 
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One  further  reason  why  the  Hertz  law  is  not  valid  is  that  for  velocities 
in  the  range  of  100  to  500  meters/sec.,  the  contact  area  diameter 
approaches  the  diameter  of  the  striking  object  , which  certainly 
violates  one  of  the  assumptions  in  the  Hertz  theory. 

At  the  writing  of  this  report,  work  is  in  progress  on  using  a 
proper  Hertz  contact  theory  for  anisotropic  materials.  However, 
for  large  deformations,  plates,  and  fracture  effects  in  the  impact 
zone,  a  more  detailed  analysis  of  the  impact  zone  is  needed  than 
is  given  in  this  report. 

Results  for  Impact  Stresses 

Solutions  to  equations  (II-2,  11-3)^  which  govern  the  central 
impact  of  anisotropic  plates,  were  found  for  impact-like  pressures 
using  an  analytical/computational  method  as  described  in  Appendix  A 
of  this  report.  In  this  model  there  are  eight  different  stresses  (see 
Section  1)  associated  with  the  membrane,  flexural  and  interlaminar 
stresses.  The  presentation  of  all  of  these  stresses  for  different 
times  after  impact  and  various  layup  angles  becomes  an  enormous  task.  Instead 
certain  key  stresses  or  stress  measures  are  presented  in  this  report 
to  give  an  overall  view  of  the  stress  picture. 

The  three  stress  measures  chosen  for  this  report  were  the  average 
membrane  stress  (tj^  +  t°^)/2  ,  the  average  flexural  stress  (tj^  +  t^g)/2 
at  the  surface  of  the  plate,  and  the  maximum  interlaminar  shear  stress 
given  by  (t^  +  •  These  stresses  are  not  necessarily  the 

maximum  stresses  at  a  point,  but  they  are  independent  of  the  orientation 
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about  an  axis  normal  to  the  plate.  The  program  also  allows  individual 
stresses  to  be  obtained  if  desired. 

The  stresses  were  calculated  in  a  quarter  plane  of  the  plate  for  a 
specific  time  after  the  initiation  of  impact  and  were  normalized  with 
respect  to  the  maximum  impact  pressure  as  calculated  in  the  above  section. 

The  data  is  presented  for  various  times  and  layup  angles  in  the  form  of 

stress  contour  plots (Fig.  14-25)  which  were  made  on  a  Calcomp  plotter 

at  the  NASA  Lewis  Research  Center,  (ref.  16) .  Superimposed  on  these  curves  are  the 

theoretical  wave  front  for  the  particular  wave  in  question  and  the 
radius  of  the  circle  which  bounds  the  impact  pressure.  The  wave 
front  calculations  were  based  on  the  work  and  represent  a 

locus  of  wave  surfaces  originating  from  the  edge  of  the  impact  circle 
for  a  given  time  after  start  of  impact. 

These  results  show  the  effect  of  the  change  of  fiber  layup 
angles  on  the  stress  distributions.  In  general,  the  maximum  stresses 
tend  to  lie  along  the  fiber  directions.  For  most  of  the  cases  treated, 
the  significant  stressed  region  revealed  by  the  calculations  is  bounded 
by  the  theoretical  wave  fronts  as  calculated  in  (ref.  1)  .  This  provides 
a  check  on  the  accuracy  of  the  approximations  made  in  the  numerical 
computation  scheme. 
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For  this  elastic  model,  the  maximum  mean  stresses  T^Ct  +  t_  ) 

■  ~  11  3  3 

occur  at  about  the  end  of  the  impact  time.  The  question  about  an 
optimum  fiber  layup  angle  is  partially  answered  by  the  data  in 
Fig.  33.  For  the  flexural  stresses,  the  optimum  layup  angle  is 
about  ±15°  showing  a  34%  lower  mean  stress  level  than  the  ±45 
case.  However,  regarding  the  interlaminar  stresses,  for  the  same 
impact  conditions,  there  seems  to  be  little  difference  in  the  maximum 
stress  level  with  layup  angle  despite  significant  changes  in  stress 
distribution  in  space  with  layup  angle. 

For  the  average  membrane  mean  stress  t°  +  t°  ,  immediately 

1133 

after  impact,  the  lower  layup  angle  plates  yield  lower  maximum 
stresses.  However,  at  later  times,  the  ±45°  layup  case  appears  to 
give  a  lower  maximum  stress  in  the  plate. 

Another  result  of  these  calculations  is  that  the  induced  stresses 
depend  on  the  impact  circle  radius- to"plate" thickness  ratio  (a/2b). 

On  the  other  hand,  the  impact  circle  depends  on  the  incoming  particle 
velocity  as  determined  by  equation  (8).  Thus,  for  each  impact  velocity, 
a  different  impact  radius/thickness  ratio  must  be  chosen  as  well  as  a 
different  contact  time.  The  integration  of  these  two  programs  has  not 

been  performed  at  this  time  but  will  be  attempted  in  the  near  future. 

Of  course,  to  evaluate  the  possibility  of  fracture  or  failure  of 
the  composite  under  impact,  the  complete  stress  matrix  at  a  point  must 
be  known,  as  well  as  the  failure  criteria,  which  will  itself  be 
anisotropic  (ref.  8) .  The  integration  of  programs  described  in  this 
section  into  a  computer  code  suitable  for  use  for  pre-design 
engineering  calculations  is  to  be  the  next  phase  of  this  report. 
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CONCLUSIONS 


The  successful  application  of  composite  materials  to  jet  engine 
fan  blades  depends,  in  part,  on  the  ability  of  these  materials  to 
retain  structural  integrity  under  transient  loadings  due  to  bird 
strikes  or  hailstone  impact  or  other  foreign  object  encounter. 

While  there  are  a  number  of  experimental  investigations  connected  with 
this  problem,  theoretical  understanding  of  impact  response  and  damage 
is  lacking.  Such  understanding  might  enable  a  reduction  in  costly 
empirical  studies  and  testing.  This  report  presents  the  first  of  a 
series  of  analytical  models  to  attempt  to  understand  impact  mechanics 
of  composite  plates. 

Using  the  method  of  stress  wave  analysis,  the  stresses  induced 
during  and  after  impact  with  a  line  load  and  central  impact  have  been 
determined.  The  model  has  been  put  into  a  computer  program  where  the 
transverse  loading  force  on  the  composite  plate  is  arbitrary.  For 
central  impact  the  results  indicate  that  the  energy  propagates  into 
the  plate  in  the  form  of  extensional  and  flexural  waves.  At  several 
impact  circle  diameters  from  the  center,  the  highest  stresses  propagate 
along  the  fiber  directions.  The  speed  of  propagation  varies  with  the 
direction  in  the  plane  of  the  plate. 

It  has  been  determined  that  for  composites  similar  to  graphite/ 
epoxy,  there  is  an  optimum  layup  angle  near  +15°  for  which  the  flexural 
stresses  are  minimized.  The  maximum  stress  in  a  plate  without  edges, 
due  to  central  transverse  impact,  occurs  at  the  end  of  the  contact  time. 
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Using  a  modified  Hertz  contact  theory,  an  estimate  of  the  contact 
times,  pressures  and  forces  for  various  impact  velocities  and  sizes 
of  ice  and  granite  spheres  has  been  calculated.  However,  the  effect 
of  the  plate  motion  has  been  neglected,  which  might  increase  the 
calculated  contact  time. 

The  next  parameters  to  study  in  this  program  are  the  effect  of 
edges  on  impact,  and  the  effect  of  plate  motion  on  the  contact  time 
and  pressures.  Also  the  validity  of  the  Hertz  theory  is  in  question 
for  the  impact  velocity  range  of  interest.  A  modified  Hertz  theory 
or  a  fluid-solid  interaction  model  should  be  developed.  Some  of  these 
effects  will  be  investigated  in  a  continuing  study  of  impact  and 
stress  waves  in  composite  fan  blades. 
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6  Center  Midplane  Displacement  vs  Time  55%  Graphite -Epoxy 
Composite,  ±45°  Layup  Angle  Impact  Load  Along  X3  Axis, 
a/b=10,0.  Comparison  of  Mindlin  and  Classical  Theories 

7  Center  Midplane  Displacement  vs  Layup  Angle  55%  Graphite 
Fiber-Epoxy  Matrix  Composite  a/b=10.0 

8  Mean  Inplane  Stress  (t°  +t°  )/2  vs  Distance  x/a  55%  Graphite 

11  33 

Fiber-Epoxy  Matrix  Composite  Plate  for  Various  Layup  Angles, 

Load  Along  X  Axis 
3 

9  Diagram  of  Plate,  Geometry  Coordinates 

a)  Wave  Surfaces  for  Graphite/Epoxy  0°  and  ±15°  Fiber  Layup 
Angle 

b)  Wave  Surfaces  for  Graphite/Epoxy  ±30°  and  ±45°  Fiber 
Layup  Angle 

11  Hertzian  Impact  Time  Versus  Velocity  for  Ice  and  Granite  Spheres 

12  Hertzian  Impact  Pressure  and  Force  Versus  Velocity 

13  Impact  Load  Distribution 

14  Average  Membrane  Stress  Contours  (tii+t33)/2  for  Graphite 

Fiber-Epoxy  Matrix,  After  Impact  for  0°  Fiber  Layup  Angle, 
Normalized  Time  t=2T 
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15 


16 


17 


18 


19 


20 


21 


22 


23 


24 


25 


Average  Membrane  Stress  Contours  (t  +t  )/2  for 

1133  o 

Graphite  Fiber-Epoxy  Matrix,  After  Impact  for  15 
Fiber  Layup  Angle,  Normalized  Time  t=2'r 

Average  Membrane  Stress  Contours  (t  +t  )/2  for 

1133  0 

Graphite  Fiber-Epoxy  Matrix,  After  Impact  for  30 
Fiber  Layup  Angle,  Normalized  Time  t=2T 


Average  Membrane  Stress  Contours  (t  +t  )/2  for 

1133  o 

Graphite  Fiber-Epoxy  Matrix,  After  Impact  for  45 
Fiber  Layup  Angle,  Normalized  Time  t=T 


Average  Membrane  Stress  Contours  (t  +t  )/2  for 

1133  Q 

Graphite  Fiber-Epoxy  Matrix,  After  Impact  for  45 
Fiber  Layup  Angle,  Normalized  Time  t=3T 

Average  Membrane  Stress  Gontours  (t  +t  )/2  for 

1133  o 

Graphite  Fiber-Epoxy  Matrix,  After  Impact  for  45 
Fiber  Layup  Angle,  Normalized  Time  t=2T 


Average  Flexural  Stress  Contours  (t  +t  )/2  for 

1133  0 

Graphite  Fiber-Epoxy  Matrix,  After  Impact  for  15 
Layup  Angle,  Thickness/Impact  Diameter  Ratio  10.0 


Average  Flexural  Stress  Contours  (t  +t  )/2  for 

1133  o 

Graphite  Fiber-Epoxy  Matrix,  After  Impact  for  30 
Layup  Angle,  Thickness/ Impact  Diameter  Ratio  10.0 


Average  Flexural  Stress  Contours  (t  +t  )/2  for 

11  33  o 

Graphite  Fiber-Epoxy  Matrix,  After  Impact  for  45 
Layup  Angle,  Thickness/Impact  Diameter  Ratio  10.0 


Average  Flexural  Stress  Contours  ft  +t  )/2  for 

11  33 

Graphite  Fiber-Epoxy  Matrix,  After  Impact  for  30^,  45° 
Layup  Angles,  Thickness/ Impact  Diameter  Ratio  =  10.0 


Maximum  Interlaminar  Shear  Contours  for  Graphite 
Fiber-Epoxy  Matrix  After  Impact  for  ±30  Fiber 
Layup  Angle,  Thickness/Impact  Radius  Ratio  =  10.0 


Maximum  Interlaminar  Shear  Contours  for  Graphite 
Fiber-Epoxy  Matrix  After  Impact  for  ±45°  Fiber 
Layup  Angle,  Thickness /Impact  Radius  Ratio  =  10. 0 
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26  Midplane  Transverse  Displacement  Distribution  in  the  Quarter 
Plane  of  the  P^ate  for  Graphite  Fiber-Epoxy  Matrix  After 
Impact  for  ±45°  Fiber  Layup  Angle,  Thickness/ Impact  Radius 
Ratio  =  1.0 

27  Midplane  Transverse  Displacement  Distribution  in  the  Quarter 
Plane  of  the  Plate  for  Graphite  Fiber-Epoxy  Matrix  After 
Impact  for  ±45°  Fiber  Layup  Angle,  Thickness /Impact  Radius 
Ratio  =1.0 

28  Average  Membrane  Stress  Distribution  in  the  Quarter  Plane 
of  the  Plate  for  Graphite  Fiber-Epoxy  Matrix  After  Impact 
for  0°  Layup  Angle 

29  Average  Membrane  Stress  Distribution  in  the  Quarter  Plane 
of  the  Plate  for  Graphite  Fiber-Epoxy  Matrix  After  Impact 
for  ±45  Layup  Angle 

30  Average  Flexural  Stress  Distribution  in  the  Quarter  Plane 
of  the  Plate  for  Graphite  Fiber-Epoxy  Matrix  After  Impact 
for  ±45°  Layup  Angle,  Thickness/Impact  Radius  Ratio  =10.0 

31  Average  Flexural  Stress  Distribution  in  the  Quarter  Plane 
of  the  Plate  for  Graphite  Fiber-Epoxy  Matrix  After  Impact 
for  ±45°  Layup  Angle,  Thickness/ Impact  Radius  Ratio  =10.0 

32  Maximum  Interlaminar  Shear  Stress  Distribution  in  the 
Quarter  Plane  of  the  Plate  for  Graphite  Fiber- Epoxy 
Matrix  After  Impact  for  ±45  Layup  Angle 

33  Flexural  Mean  Stress  History,  t  +t  versus  Time  After 

11  33 

Impact,  at  Origin  for  Graphite-Fiber  Epoxy 


34  Interlaminar  Flexural  Stress  History  /t^  +t^ 

12  32 

Versus  Time  After  Impact  for  Graphite-Epoxy  r/a=0.6 


Note : 


T  =  Time/a,  (ysec/mm] 

for 

Figs. 

T  =  Time/[b/(C  /p)l/2)] 

for 

Figs  . 
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APPENDIX  A:  THEORETICAL  ANALYSIS  OF  IMPACT  IN  COMPOSITE  PLATES 


Introduction 

In  a  previous  paper,  (ref.  I)  and  in  section  I  of  this  report 
the  Author  examined  the  propagation  of  wave  surfaces  in  composite 
plates,  and  the  response  to  a  line  impact  load  respectively.  In 
section  II  of  this  report  results  for  the  two-dimensional  plate 
response  to  a  distributed  impact  load  were  presented.  This  analysis 
makes  use  of  a  computational  tool  called  the  "Fast  Fourier  Transform," 
which  permits  the  calculation  of  the  inverse  of  Fourier  transforms  on 
the  digital  computer.  The  application  of  this  technique  to  the  calculation 
of  impact  induced  stresses  is  renewed  in  this  section. 

The  mathematical  model  treats  the  laminated  plate  as  an  equivalent 
anisotropic  material  using  a  program  developed  by  Chamis  (ref.  9) .  A  modi¬ 
fication  of  Mindlin's  theory  of  crystal  plates  is  used,  which  results 
in  five  two-dimensional  stress  waves.  Two  of  these  waves  describe  the 
average  or  membrane  stresses,  while  three  other  waves  are  associated 
with  the  flexural  motion.  The  two  former  waves  are  non-dispersive , 
while  the  flexural  waves  exhibit  strong  dispersion  at  the  low  frequencies. 
In  the  Mindlin  plate  theory,  (ref.  7),  the  impact  transverse  surface  force 
enters  the  problem  directly  through  the  differential  equations.  To  solve 
this  initial  value  problem,  a  Laplace  transform  is  taken  on  the  time 
mriable,  and  a  double  Fourier  transform  is  taken  on  the  two  space 
variables.  A  solution  of  this  triple  transformed  problem  is  obtained 
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in  the  transform  space.  Finally,  the  solution  is  completed  by  an 
analytical  inversion  of  the  double  Fourier  transform  using  the  fast 
Fourier  transform  algorithm. 

The  displacement  variables  used  in  the  theory  are  described  as 

follows:  u°  and  u*^  represent  the  midplane  displacements  in  the  plane 

of  the  plate,  (see  Fig.  9);  u°  represents  the  transverse  midplane 

2 

displacement  and  u^  are  a  measure  of  the  rotation  of  a  line 

1  3 

normal  to  the  midplane.  The  equations,  which  were  derived  in  Reference  1, 
are  listed  below,  where  is  a  transverse  tension  force  on  the  plate 

surface  and  the  constants  are  the  equivalent  elastic  stiffnesses 

for  the  composite  plate. 
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where 


^11  “  ^ir  ^12^^22 


S3  =  S3-  S2/S2 


*"13  “  *"13  *"12*"32/*"22 


Se  =  ^  Se 


K  =  1x2/12 

Analytical  Part  of  Solution:  Midplane  Motion 


Solutions  for  the  transforms  of  u*^ 

1 

easily  found.  Tlie  Laulace  variable 


and  u*^  are 
3 


The  Laplace  variable  is  denoted  by  s  and  the 


Fourier  variables  are  =  k  cos  a,  k^  =  k  sin  a  .  The  vector  formed 
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from  (k  ,  k  )  =  k  represents  a  Fourier  wave  number  vector  corresponding 
to  the  harmonic  frequency  u  =  -is  .  The  resulting  expressions  become. 


T[L[u«l] 

1 

=  u  = 

1 

i  q 

{CA22  + 

s^/k^) a^ 

-  Ai2a2 

sin^a} 

cos 

a/kA 

T[L[u»]] 

3 

& 

=  u  = 

3 

A 

i  q 

{(A^l  - 

s2/k2)a2 

-  N2^1 

cos^a} 

sin 

a/kA 

CA-3) 

A  = 

(A^l 

+  s^/k^) 

(Aj2  •  S 

Vk2)  - 

A^  cos 
12 

sin^a 

where , 

A 

C,,  sin^a  +  C^.cos^a, 


3-1  -  C12/2C22  >  ^2  “ 

(A  bar  indicates  a  Laplace  transform,  and  q  indicates  a  double  Fourier 
transform.)  These  solutions  have  the  form. 

0  =  i  q  P(s2)/kA(s2) 

where  P(s^),  A(s^) ,  are  polynomials  in  s^  .  ^(s^)  has  four  zeroes 

in  the  complex  s  plane  for  each  1^.  These  roots  have  the  form, 

s  =  ±  ivj^  k  ,  ±  iv2  k 


cos^a  +  CggSin^a  ,  A 


22 


^12  "  ^13  ^55 
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where  v  and  v  are  the  plane  wave  velocities  corresponding  to  the 
1  2 

wave  normal  k  .  Values  o£  v  and  v  versus  a  were  reported  in 
-  1  2 

Reference  1.  The  Laplace  inversion  of  (A-3)  can  then  be  done  by  use  of 
the  convolution  theorem; 


U(k,a,t)  =  ^ 


q(k,a,t-T)  G(T)dT 


(A-4) 


GCt)  =  - 


k  P(-k2v2)  sin  k  v^t  k  P(-k2v|)  sin  k  v^t 


(v^-v^ 


V, 


(v^-V^) 


Analytical  Part :  Flexural  Motion 

By  a  similar  procedure,  one  can  solve  for  u°  ,  u|  ,  u^  ,  which 
have  the  form 


W  = 


.  R(£) - 

(s2,k2) 


(A-5) 


However  in  this  case  the  roots  of  =0  are  not  proportional  to 
k  ,  This  means  that  the  velocity  of  the  waves  depends  on  the  wavelength. 

It  is  known  from  the  dispersion  relations  for  these  plates  (ref.  7) ,  that  for 
k  real,  there  will  be  six  pure  imaginary  roots  of  A^ (s^)  =  0; 


s  =  iiw^Ck),  ±ia)2(k),  iiwgCk) 

The  Laplace  inversion  of  (A-5)  again  makes  use  of  the  convolution 
theorem  and  the  residue  calculus  to  invert  1/A(s2).  Thus  we  obtain 
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W(Jc,t)  =  I  q(k,t-T)  H(T)dT  (A-6) 

■'  0 

where 

Rfo),)  sin  (io,t  R(ti)  )  sin  to,t 

H[t)  =  — - ! -  "  — ^ - ^ 

(oj2-a)p  (a)2-a)2)(j^  (to2-a)2)a)^ 


RCcOg)  sin  ojgt 
(0)2-0)^ 


Numerical  Inversion: 

The  inversion  of  the  Fourier  transforms  involves  integrals  of 
the  form 


U(x^  jXg) 


4tt^ 


U(k^,k2)  e 


-")s-«dk 


dk. 


CA-7) 


If  significant  changes  in  UCx)  occur  over  distances  greater  than  X  , 
then  the  largest  wave  ntimber  of  interest  will  be 

K  =  27r/X  . 

Thus  we  may  be  satisfied  with  an  approximation  to  UCx)  of  the  form 


U(^)  = 


4'n'^ 


U(k)  e"^^'^dk,dk. 

1  3 


Shifting  coordinates ,  this  becomes 


fit?)  = 


^iCKx^+Kxg) 

4ir^ 


2K 


?)(k 


^-K,kg-K)  e  ^v'^dk^dkg 
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This  latter  integral  may  further  be  approximated  by  the  following  sum. 


U  = 


^iCKx^+KXg)  (j3_K3  3-i[k^(J)xi+k^(J)x^] 

7r2N2  I,J=1  ^  ^ 


where 


k^d) 


iL 

N  2 

. 


(I-l) 


k  (J) 

3 


2K 

N 


I  *  CJ-i) 


or 


it2n2 


I ,J=1  (A- 8) 


When  X  and  x  are  continuous  variables,  a  summation  of  the  type 
1  3 

above  must  be  performed  for  each  point  (x^,X2),  which  makes  the  cal¬ 
culation  of  the  sums  impractical  for  a  large  number  of  grid  points 

fx  X  1  However,  if  x  and  x  take  on  certain  descrete  values,  there 
''1*3  1  3 

is  an  algorithm  (ref.  8) ,  which  makes  the  calculation  of  these  sums  feasible  for 
a  large  number  of  grid  points.  This  algorithm,  known  as  the  fast  Fourier 
transform,"  takes  a  sampling  matrix  of  the  transform,  say  A(kj,kj),  and 
returns  the  following  sampling  matrix  of  the  original  transformed  function. 


T(L,M)  = 


E 

J=1 


N 


E 

J=1 


A(k, 


,kj) 


,-2,i 


CJ-1)(M-1) 

N 


■] 


(A-9) 


This  operation  is  known  as  a  descrete  finite  Fourier  transform.  To 
put  the  above  expression  into  this  form  we  choose  for  the  discrete 


values  of  x.  and  x  the  numbers 
^  3 
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Xj(L)  =  I  (L-1)  ,  x^CM)  =  I  CM-1) 


(A-IO) 


Numerical  Results 

Several  checks  were  made  of  the  fast  Fourier  computer  routine 
(ref.  8)  used  in  this  paper.  First,  known  functions  were  transformed 
analytically  and  inverted  numerically.  These  tests  revealed  that 
one  must  work  with  continuous  fimctions  if  spurious  oscillations 
near  points  of  discontinuity  are  to  be  avoided. 

Second,  known  one-dimensional  solutions  were  checked  using  the 
numerical  transform  method,  such  as  the  impact  of  a  string  and  the 
impact  of  a  classical  beam  (see  Section  I).  All  these  checks 
revealed  very  close  agreement  between  the  numerical  transform  output 
and  the  known  analytical  result. 

The  impact  pressure  distribution  used  was  the  following  (see 
Figure  13) 


q,  =  -P„  (1  -  2  (^)  +  (^)  )  sin 


Trt 


for  r<a,  x^  +  x|)  and  t<T 

i  o  0 


(A-11) 


q„  =  0,  for  r>a  or  t>T 
^  o 

The  Hertz  contact  pressure  based  on  static  isotropic  elasticity  has  the 
form 
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CA-12) 


The  principal  difference  between  (A-11)  and  (A-12)  lies  in  the  infinite 
slope  in  (A-12)  at  r=a  .  The  form  of  the  impact  pressure  (A-11)  was 
chosen  to  avoid  such  infinite  slope.  A  distribution  with  continuous 
first  and  second  derivatives  (as  the  function  (A-11)  exhibits)  was  dictated 
by  a  desire  to  have  all  stresses  continuous  (i.e.  to  avoid  shocks)  and 
thus  avoid  spurious  oscillations  in  the  numerical  inversion.  This  require¬ 
ment  results  from  the  fact  that  in  the  Mindlin  theory  the  average  mid¬ 
plane  stresses,  having  as  wave  sources  terms  proportional  to  Vq2 ,  have 
the  form 


t  „  [  kS  (kT)sin  k  v(t-T)dT  (a,3  =  1,3) 

a3  ^ 

where  v  denotes  either  the  first  or  second  wave  speed.  Thus  when 
q^  has  the  form  of  a  short  duration  impact  i.e. 

=  Q(k)6(T) 

the  Fourier  transforms  of  the  stresses  are  proportional  to 

t  „  k2Qrk)  sin  k  vt 
a3  - r - 
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or  for  small  times. 


9r2 

Thus  for  the  stresses  to  be  continuous  at  early  times  after  the  wave 
arrival,  the  spacial  part  of  must  have  continuous  second  derivatives, 

which  led  to  the  choice  of  (A-11) .This  conclusion  was  also  reached  in  numeri 

cal  tests  of  the  fast  Fourier  program  when  non-smooth  load  distributions 

were  used.  (See  Section  II  of  this  report). 

Using  a  "fast  Fourier"  computer  routine  written  in  Fortran  IV, 

the  induced  impact  stresses  were  calculated  on  both  an  IBM  360-91,  and 

IBM  7094  computers .  The  grid  used  was  32  x  32  or  1024  points  in  the 

-  Xy  plane.  Execution  time  on  the  360-91  was  of  the  order  of  12 

seconds,  and  22  seconds  on  the  7094.  A  denser  grid  of  64  x  64  was 

also  tried  with  a  running  time  of  less  than  a  minute  on  the  IBM  360-91. 

The  output  data  consisted  of  a  matrix  (32  x  32)  of  stress  values 

for  the  quarter  plane  of  the  plate.  Interpolation,  contour  plotting, 

and  three-dimensional  plot  routines,  developed  at  the  NASA  Lewis  Research 

Laboratory  for  use  with  a  "Calj^  Comp"  plotter,  were  used  to  obtain  stress 

contours  and  three-dimensional  plots  as  shown  in  Figure  2  of  Section  II. 

The  significant  stress  levels  all  lie  within  the  surface  bounded 

by  the  theoretical  wave  surface.  In  Figures  14,  17  the 

average  or  membrane  mean  stress  contours  i(t'^  +  t°  )  for  graphite 

2  11  33  ' 

fiber/epoxy  matrix  laminate  plates  are  shown  for  layup  angles  of 
0°,  ±45°. 
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The  stresses  shown  correspond  to  the  faster  wave  speed  (quasi  com- 
pressional  wave)  which  is  anisotropic,  as  is  seen  from  the  wave  surfaces 
in  Figure  10  .  The  stresses  associated  with  the  second  wave  speed 
(or  quasi  shear  wave)  were  much  smaller  than  those  in  the  faster  wave 
and  are  not  presented. 

The  flexural  or  bending  motion  has  three  waves  associated  with  it. 
The  largest  stresses  however  were  found  in  the  lowest  flexural  wave 
which  travels  at  an  isotropic  wave  speed  given  by 


V3  =  [Cgg 


(k  =  77^/12  ,  is  Mindlin's  correction  factor  (I'ef-  3).  Stress  contours  for 
the  mean  flexural  stress  ^  +  ^33  in  this  wave  are  shown  in  Figures  20- 

25  for  graphite  fiber/epoxy  matrix  laminate  plates  under  the 
transverse  impact  pressure  (A-11) .  Note  that  the  wave 
front  is  circular  since  V3  is  isotropic  for  laminate  plates.  Stresses 
in  the  second  and  third  flexural  waves  were  found  to  be  small.  Three- 
dimensional  computer  plots  are  shown  in  Figures  26-32. 

The  maximum  stress  levels  were  found  to  occur  immediately  after 


the  end  of  impact  and  appeared  to  propagate  along  the  fiber  directions, 
given  by  the  layup  angles. 


APPENDIX  B:  LIST  OF  SYMBOLS 


a 


a. 

1 

A.  . 
II 

b 

C.  . 

D(I) 

e.  . 

ij 

E 


F 


M 


P 

n 

S 


t 


t.  . 

II 


u 

u 


0 

2 

1 

1 


u; 


half  width  of  impact  contact  line  (one  dimensional)  or  radius 
of  impact  contact  circle  (two  dimensional) 

ratios  of  elastic  constants 

acoustic  tensor  components 

half  thickness  of  plate 

anisotropic  elastic  constants 

column  matrix  of  discrete  Fourier  Transform 

elastic  strain  tensor 

Young's  Modulus 

impact  force 

wave  number 

Hertz  contact  law  constant 

Laplace  transform 

mass  of  impacting  object 

wave  normal 

Legendre  polynomials 

impact  loading  function 

position  vector 

Laplace  transform  variable 

time 

stress  tensor 
displacement  vector 

inplane  displacements 
transverse  plate  displacement 
flexural  rotation  displacements 
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higher  order  plate  displacements 

V  wave  velocity 

V  velocity  of  impacting  object 
cartesian  coordinates 


a  angle  of  one-dimensional  wave^normal 

a  relative  approach  of  impacting  object  and  plate 

A,A^  characteristic  acoustic  determinant 

fiber  layup  angle 

K  Mindlin  correction  factor,  or  wave  number  (see  the  text) 

A  wave  length 

V  Poisson's  ratio 

p  density 

frequency 

T  impact  time  or  time  parameter 

5  one-dimensional  wave  normal  direction 
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APPENDIX  C:  IIATERIAL  PROPERTIES 


All  the  calculations  in  this  report  are  for  the  composite  material 
consisting  of  graphite  fibers  in  an  epoxy  matrix. 

The  values  of  the  elastic  constants  for  graphite/epoxy ,  for  various 
ply  layup  angles,  were  obtained  from  an  analysis  by  Chamis  (ref.  9). 

The  assumed  properties  of  the  graphite  fibers  and  epoxy  used  by  Chamis 
are  as  follows; 

Epoxy  Young's  Modulus  E  =  0.57  lO^psi 

Poisson's  Ratio  v  =  0.36 

Graphite  Fiber  (Thomel  50) 

"1"  Axis  along  the  fiber 

E  =50  lO^osi 
1 1 

E  =  E  =1.0  lO^psi 

22  33 

V  =  V  =  0.20 

12  13 

V  =0.25 

23 

G  =  1.3  10^ 

12 

G  =  0.7  10^ 

23 

The  values  of  the  elastic  constants  for  the  composite  are  given  in 
the  following  table. 
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FIGURE  I  •  DIAGRAM  OF  COMPOSITE  PLATE 
GEOMETRY  AND  NOTATION 


WAVE  NUMBER  kb 


FIGURE  2 -DISPERSION  RELATION  FOR  FLEXURAL  MOTION  (ub  vs  kb 
FOR  55%  GRAPHITE- EPOXY  ±  45®  LAYUP  ANGLE 
FOR  WAVE  NORMAL  a  =  0  ALONG  X,  AXIS 


Q 

kJ 

N 


O' 


u 

a> 

tA 

>3^ 


£ 

£ 


z 

kJ 

5 

kJ 

y 

0. 

CO 


q: 

< 

_j 

ID 

CD 

Z 

< 

q: 


< 

H 

Z 

§ 

oc 

o 

k. 

zg 

^  X 

LlI. 

k-X 

il 

X 


11 

o 


ro 

kJ 

X 

ID 

O 

k. 


O 

z 

X 

f~ 

(A 

kJ 

t 

Z 

k. 

Z 

Z 

< 


CO 

kJ 

I 

H 

Z 

kJ 

yoD 

<  w 

-J  . 
Q_  D 
CO  O 
QkJ 

U_  CD 
OZ 

z9 

o< 

n:  o 

^kJ 

3| 

kJ 


0. 

:e 

o 

CD 


O 

I- 

kJ 

3 

Q 


o 


o 


d 


q°d/9®3  gn  lN3^N3^V^dSla  a3ZnVlNdON 


-0.2 

FIGURE 


FIGURE  5- MIDPLANE  DISPLACEMENT  WAVES, 557o  GRAPHITE- EPOXY  COMPOSITE,  ±45®  LAYUP 
ANGLE  a/b  =  10.0 ,  IMPACT  LOAD  ALONG  X3  AXIS 
COMPARISON  OF  MINDLIN  AND  CLASSICAL  THEORIES 


NORMALIZED  DISPLACEMENT  UjCse^Pob 


2  4  6  8  10 


FIGURE  6  •  CENTER  MIDPLANE  DISPLACEMENT  vs  TIME  55% 

GRAPHITE  -  EPOXY  COMPOSITE  ,  ±45°  LAYUP  ANGLE 
IMPACT  LOAD  ALONG  X,  AXIS,  a/b=  10,  COMPARISON 
OF  MINDLIN  AND  CLASSICAL  THEORIES 
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FIGURE  7- CENTER  MIDPLANE  DISPLACEMENT  vs  LAYUP  ANGLE  55% 
GRAPHITE  FIBER -EPOXY  MATRIX  COMPOSITE  a/b  =  10.0 


NORMALIZED  MEAN  MEMBRANE  STRESS 
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=  0.2  ^sec/mm 
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FIGURE  8- MEAN  INPLANE  STRESS  (tj;+t03)/2  vs  DISTANCE  x/a  55%  GRAPHITE 
FIBER -EPOXY  MATRIX  COMPOSITE  PLATE  FOR  VARIOUS  LAYUP 
ANGLES,  LOAD  ALONG  Xg  AXIS 
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FIGURE  9  CENTRAL  IMPACT  GEOMETRY 


Distance  from 
impact  point, 
mm 


Distance  from 
impact  point, 
mm 


FIGURE  10.  a)  WAVE  SURFACES  FOR  GRAPHITE/EPOXY  0°  and  ±15°  FIBER 
LAYUP  ANGLE 


t  =  10-®  sec. 


FIGURE.  10.  b)  WAVE  SURFACES  FOR  GRAPHITE/EPOXY  ±30°  and  ±45°  FIBER 
LAYUP  ANGLE 


t  =  10”^  sec. 
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IMPACT  VELOCITY,  METERS/sec 

HERTZIAN  IMPACT  CONTACT  TIME  FOR  GRANITE  AND  ICE  SPHERES 

ON  GRAPHITE/ EPOXY  COMPOSITE 


FinURE  11 


MAXIMUM  IMPACT  PRESSURE  10^  psi 


PRESSURE,  10®  NEWT./m 


MAXIMUM  IMPACT  FORCE  10^  NEWTONS 


VELOCITY  ,  METERS  /sec 


FiGURK  12b  MAXIMUM  HERTZIAN  IMPACT  FORCE 


IMPACT  FORCE  ,  lO^lbf 


MEMBRANE  STRESS  f  (t,,  +  tja) 
GRAPHITE  FIBER/ EPOXY  MATRIX 
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FIGURE  15 
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MEMBRANE  STRESS  j  (tn  +  taj) 
GRAPHITE  FIBER/ EPOXY  MATRIX 


FIGURE  17 
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FIGURE  19 


FIGURE  20 
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FIGURE  21 
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FLEXURAL  STRESS  fltu+tsa) 
GRAPHITE  FIBER/EPOXY  MATRIX 


FIGURE  23 


INTERLAMINAR  SHEAR  STRESS  ^A■(tf2+i^ 
GRAPHITE  FIBER/ EPOXY  MATRIX 


FIGURE  24 


FIGURE  25 


TRANSVERSE  DISPLACEMEN 
GRAPHITE  FIBER/ EPOXY  MA‘ 


T.'P. 

O  UJ 

C0  5 


(O 

LiJ 


(T. 


lO 


II 

u.  I- 

3  5^’ 

5-^0= 

o 


pUJ 

io 

“2  <5 

+1  •  H 
.  .CO, 
iijan 

d§§ 

2x5 

bi 

3li 


(M 

=) 


FIGURE  27 


MEMBRANE  STRESS  i(t,|+  tjj) 
GRAPHITE  FIBER/ EPOXY  MATRIX 


FIGURE  28 


Mean  In-Plane  Stress  Profile  After  Impact 


FIGURE  29 


FLEXURAL  STRESS  ?  (t„  +  t33) 
GRAPHITE  FIBER/EPOXY  MATRIX 
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FIGURE  30 
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FLEXURAL  STRESS  5(t„+ W 
GRAPHITE  FIBER/EPOXY  MATRIX 


FIGURE  32 
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NORMALIZED  CONTACT  TIME  TV3/b  =  I.O  (V3=I.I8  mm/^isec) 
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FLEXURAL  STRESS  tigj  vs  TIME  FOR  a/b=I.O  , 

NORMALIZED  CONTACT  TIME  rV3/b  =  1.0  (V3=  1.18  mm/^sec 


